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The carbon wall oxidation technique coupled with a CFD technique was employed to 
study the flow in the expanding crack channel caused by the oxidation of the channel 
carbon walls. The recessing 3D surface morphing procedure was developed and tested 
in comparison with the arcjet experimental results. The multi-block structured adaptive 
meshing was used to model the computational domain changes due to the wall recession. 

Wall regression rates for a reinforced carbon-carbon (RCC) samples, that were tested in 
a high enthalpy arcjet environment, were computationally obtained and used to assess 
the channel expansion. The test geometry and flow conditions render the flow regime 
as the transitional to continuum, therefore Navier-Stokes gas dynamic approach with the 
temperature jump and velocity slip correction to the boundary conditions was used. The 
modeled mechanism for wall material loss was atomic oxygen reaction with bare carbon. 

The predicted channel growth was found to agree with arcjet observations. Local gas flow 
field results were found to affect the oxidation rate in a manner that cannot be predicted by 
previous mass loss correlations. The method holds promise for future modeling of materials 
gas-dynamic interactions for hypersonic flight. 

I. Introduction 

The reinforced carbon-carbon (RCC) nose cup and wing leading edge of the Space Shuttle, studied 
experimentally and numerically in this research, are vulnerable to cracks that can be formed through micro- 
meteoroid and orbital debris impact. * 1 During the Earth atmosphere re-entry, the free stream species, in- 
cluding reactive atomic 0, may penetrate through the bare carbon-crack channel and oxidize the channel 
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walls. The resulting growth of the crack can pose severe structural damage to the wing and the nose portion 
of the vehicle. The dimensions of the cracks are sufficiently small while the flow presiding the channel has 
larger characteristic sizes. This is why for in-flight re-entry conditions, as well as ground based testing, the 
multi-scale flow solutions should be obtained. 

The crack expansion process in RCC samples, modeling a portion of the Space Shuttle wing, has been 
studied in ground-based arcjet tests 2 and the high velocity meteoroid impact holes subjected to the arcjet 
flow conditions were presented in. 3 The purpose of this work is to demonstrate a physical model that may 
be employed in numerical simulations of microflows through a carbon channel and to predict the extent 
of the channel expansion in the geometrically complex configurations including 3D Space Shuttle leading 
edge fragment and axisymmetric stagnation region hole which includes a sub-channel causing delamination. 
Atomic oxygen impacting a bare carbon wall was assumed to chemically react with the channel walls to 
produce CO gas that then enters the flowfield. The probability of the reaction is the only free parameter in 
the model, with a value of 0.9 used in all simulations. 4 

II. Experiment Conditions 

The high energy flow is generated in the arcjet facility through an expansion of the heated and pressurized 
gas through a supersonic nozzle. Figure 1 presents the stagnation point hole experimental setup, which is 
referred to as Case 1151 in this paper and Figure 2 presents the schematics of the experimental facility along 
with the major flow and device parameters as well as sub-zones for NATA and CFD (Gasp) 5 solutions for 
the case referred to as Case 2033. The specimen presented on the right hand side of the Fig. 2 resembles 
a portion of the Space Shuttle leading edge with a tiny channel machined through the surface which is 
exposed to the flow. The initial (left hand side) and the post-test (right hand side) configurations of the 
specimen Case 2033 are presented in Fig. 3. It is clearly visible in the figure that the initial micro-channel 
has significantly expanded, when subjected to the arcjet flow conditions similar to that during the vehicle 
Earth atmosphere re-entry. Similar observation is evident for the case 1151. If happened during the flight, 
the damage of such a magnitude may severely compromise the wing structure by exposing its interior to 
the high energy re-entry flow. The goal of this research is therefore to develop and test a computational 
technique which can predict the microchannel expansion when subjected to the chemically active re-entry 
flows. 


III. The Computational Technique 

In order to model the channel expansion, a CFD solver should be enabled with the specific boundary 
conditions modeling the chemical processes at the walls of the channel. Because of its robustness, when used 
for the solutions of the high and low speed compressible flows, the Gasp 5 N-S based solver was used in this 
research. An additional type of boundary conditions previously not available in GASP was implemented 
based on the below explained formulation. 
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An important feature of the crack flow modeling is the interaction between the gas and the channel wall, 
which leads to mass, momentum and energy exchange. It was found that the chemical reactions among gas 
species (such as O + N 2 — > NO + N) are rare in the channel flow since the total reaction energy of the 
colliding molecules for a gas temperature of 4000 K is usually less than the chemical activation energy. As 
is discussed in Ref., 6 atomic oxygen plays an important role in carbon mass surface loss. Reference 4 gives 
a probability close to 0.9 for the gas surface reaction, 


C (s) + O CO 


(1) 


where, the temperature of gaseous atomic O is about 1500 K, C( s ) is a carbon atom at the wall surface, 
and the other species are gaseous. We can estimate the mass loss due the above gas-surface reaction in a 
simple calculation. As the gas travels through the channel, it interacts with the carbon wall. The number 
flux normal to the channel wall, N, is given as, 7 


N = 



(2) 


where no is the atomic O number density and /? is the reciprocal of the molecular most probable thermal 
velocity. The mass loss, surface regression rate, S, is related to N by, 

N 

S=—P (3) 

nc 

where nc is the carbon number density (a mass density value of 1.575 gm/cm 3 was assumed) and P is the 
reaction probability for the C( s ) + O reaction. Using typical values in the boundary layer region of the 
channel flowfield of no = 3.0 x 10 21 m -3 , T = 2400 K, and a reaction probability of 0.9, we obtain an 
approximate surface regression rate of 0.02 mm/s. This wall regression rate will be shown to be typical of 
the values measured in the arcjet testing. 

The typical magnitude of the wall regression rate, when compared with the channel dimensions (from 1 
to 50 mm), suggests that the material response and the gas dynamic flows may be loosely coupled, thereby 
simplifying the crack wall damage modeling. Since it takes the gas less than 10 -3 sec to travel through the 
channel, it is reasonable to assume that the flow through the channel reaches a steady state at each of the 
various, time-dependent channel geometries. 

Modified boundary conditions which account for the processes given by Eq. (1) were developed for the 
GASP commercial CFD software which solves N-S equations. The numerical solution of the Navier-Stokes 
equations for viscous flow was obtained by use of the finite volume discretization of the computational 
domain done with the Gridpro 8 meshing software. The temperature jump and velocity slip corrections to 
the boundary conditions were used to account for the rarefaction effects inside the channel. 
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IV. Computational Domain and expanding surface morphing 

A. Topology based meshing 

In general, block structured meshing is a well defined notion - meaning that the entire computational domain 
is divided into blocks and each of the blocks has a structured mesh. Since topology is a more abstract notion 
let us define it. 

In the geometrical sense, topology is defined as an organization of space. For example one can think of 
geometrical figures as topology elements, but not all geometrical figures represent distinct topology types. 
If we take for example a torus with some size as a topology example and stretch or bend it, or do whatever 
else we want to do with its shape and characteristic sizes, the resulting geometrical figure will not represent 
a new topology as long as our movements do not make the torus surfaces intersect each other. In this sense, 
a torus and say a picture frame are topologically equivalent. The concept behind Gridpro 8 meshing rests on 
the idea that a topology remains unchanged even though the shape of the corresponding geometrical entities 
(blocks of mesh in this case) change to some extent as explained above. Topology usage for meshing has 
some specifics. Let us cover some of them. 

In order to mesh the volume regions, one must cover them with wireframe topologies that approximately 
follow their surface boundaries. The majority of Gridpro 8 operations deal with creating and modifying 
topologies. We note that such topologies consist of a wireframe and surface assignments. The wireframe 
consists of topology nodes and edges. The assignment process is a way for telling the Gridpro 8 software to 
which particular boundary one or another piece of topology belongs to. 

The next issue is how accurately the user should position the topology elements (nodes and edges). 
Gridpro 8 can work with loosely positioned topologies. In other words, it can grid a topology which is 
defined with nodes and edge positions that only roughly conform to the desired (corresponding) geometry. 
A wireframe topology does not have to exactly fit the boundaries for the region to be meshed. However, 
the topology itself (the way the topology elements are assembled and assigned to the surfaces) must not 
have errors and if an error in the topology is made the meshing will not proceed. If a topology element 
(a node or an edge) is positioned too far from a surface, this will not create a fatal error but the meshing 
process will take longer than it otherwise should. Thus, as long as the topology is positioned with some 
care, the surfaces and volumes will be meshed with relative ease. There is no a priory way to tell how close 
the topology elements should be placed with respect to the corresponding surfaces. The trade off is that the 
closer a topology resembles the actual geometry, the faster the meshing procedure will converge. The final 
mesh quality does not depend on how close the topology elements are located to the corresponding surfaces, 
but again, errors in the topology, or a topology which does not properly reflect the underlying computational 
domain, may result in a fatal problem during the meshing procedure or a mesh of poor quality. 

Gridpro 8 grids flow volumes in both two and three dimensions and therefore topologies can be two and 
three dimensional in Gridpro. The 1151 case is solved in two dimensions while 2033 is a full 3D case. However 
some initial modeling of the 3D problems as 2D may help building a better topology and the resulting mesh. 
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In two dimensions, gridded blocks are represented by convex quadrilaterals, or quads and in three dimensions, 
griclded blocks are represented by convex hexahedrons, or hexes. 

After this short introduction to the topology as it is used in Gridpro, let us consider the 1151 case and 
decompose its computational domain into topology quads and then mesh the domain. Because surfaces can 
be complex, they may require complex topologies and unfortunately no unique topology exist for a case. 
Furthermore different topologies will produce meshes of different quality. Let us now start with the simplest 
topology possible for the case 1151 and see the results. Figure 4 presents such a topology consisting of 18 
topology elements and 14 topology nodes. As can be seen the topology presented in Fig. 4 roughly mimics 
the domain boundaries and physical surfaces presented in Fig. 5, however two observations should be done. 
The first one is that the computational domain represented by the 2D surfaces in Fig. 5 exactly describes the 
specimen and the outer domain boundaries while the topology presented in Fig. 4 has a similar shape but 
its edges and nodes do not necessarily lay on the surfaces. The second observation is that even this simple 
topology contains some additional edges and nodes which are not required to define the domain boundaries 
and body surfaces. These edges and nodes are created keeping in mind the extended structure of the mesh 
blocks and are in fact required. The topology nodes are assigned to the corresponding surfaces, which in 
this case is rather obvious and one proceeds with the meshing. Recall that assignment of a node to an edge 
is used to connect the wireframe topology to the true physical surfaces of the computational domain. 

Figure 6 presents the resulting mesh block structure obtained after a few hundred steps of the meshing 
procedure solving the variational equations as defined by the underlying algorithm. 8 As can be seen in the 
figure the resulting block structure is topologically equivalent to the wireframe structure presented in Fig. 4, 
however the outer boundaries of the mesh blocks now exactly matching the boundaries of the 1151 specimen 
and the computational domain. The inner boundaries of the mesh blocks have also been modified to create 
an optimum mesh structure based on the provided topology. 

Let us take a look at the mesh itself. The mesh density is kept relatively coarse for this illustrative 
example to reduce clutter. Figure 7 presents the entire mesh for the case 1151 with an enlargement of the 
channel area presented in Fig. 8 and the corner of the channel entrance shown in Fig. 9. As can be seen 
in Fig. 7, the overall mesh quality is reasonably good with acceptable meshing inside the channel. Let us 
now look into some details. It can be seen in Fig. 9 that the mesh in this area is actually of a very poor 
quality. None of the quality criteria such as orthogonality or skewness can be satisfied with such a mesh. In 
fact if one attempts to use such a mesh for a CFD solution, the solver would probably fail or produce some 
unphysical results at best. 

Let us see what can be done to improve the mesh quality. We start with analyzing what causes the 
problem which arises from the so called ’’singularity” at the corner of the channel. There are several types 
of singularities defined in Gridpro (which should not be confused with the mathematical definition of a 
singularity). In this case, the particular type of singularity arises because four topology edges and only two 
surface edges are converging at the corners of the channel. 

Let us try avoiding the problem. To avoid it we define a new 1151 topology shown in Figure 10 with 


5 of 41 


an enlargement of the channel area presented in Fig. 11. Figure 12 presents the new surfaces defining 
the computational domain. There are several differences between the previous and the new topology and 
corresponding surfaces involved in the domain modeling. The first difference is the addition of the so called 
artificial surfaces which can be observed in Fig. 12. These surfaces do not represent real boundaries of the 
computational domain. Their sole purpose is to resolve the singularity issue at the corners of the specimen. 
Their presence however must be reflected in the topology structure for the effect to take place. As can be 
seen in Figs. 10 and 11, several new topology blocks have been added corresponding to the new domain 
model. The topology nodes must now be partly reassigned to the new surfaces which came into picture. 
Some of the nodes must now have double or tripple assignments. For example the topology node located at 
the corner of the channel entrance (see Fig. 11) is now assigned to three surfaces: to the channel wall, the 
specimen outer surface, and to the artificial surface located at the entrance of the channel. 

Let us now look at the resulting mesh. The block structure of the new 1151 mesh is presented in Fig. 13. 
It is immediately clear from the observation of the block structure that the meshing problem at the corner of 
the channel entrance is now resolved. The angles at the corners of the blocks are close to 90 degrees, which 
is a good indication of the mesh orthogonality in the problematic areas. The boundaries of the blocks are 
also adapted to produce a better mesh based on the new version of the topology. Let us take a look at the 
mesh itself. Figure 14 presents the entire mesh for the the domain while Fig. 15 presents the mesh inside the 
channel with an enlargement of the channel corner presented in Fig. 16. As can be seen in the figures, and, 
in particular in Fig. 16, the mesh quality is improved dramatically and this new mesh can definitely be used 
for the CFD computations. One additional implication of using the artificial surfaces is the improved ability 
to manipulate the local mesh density. In the 1151 case we want the mesh to be dense inside the channel to 
resolve the area of interest and less so in the outer portion of the domain to save computational cost. With 
the artificial surfaces in place such a mesh can be produced. The artificial surfaces bound the interior of 
the channel preventing the mesh optimizer from moving the mesh points located inside the channel into the 
areas in front and behind the channel. 

With the acceptable mesh having been created we conclude the introduction of the topology concepts. 
Our resulting mesh may still need additional refinement, however. Issues such as the shock wave resolution 
in front of the specimen and optimization of the mesh density need to be solved before the actual CFD 
calculations may be done. All of these issues are addressed and will be commented in result section of this 
article. Cases 1151 with delamination and 2033 solved here are significantly more involved, but the presented 
concepts hold. Let us now introduce the additional subjects of 3D surface definitions and surface morphing, 
relevant to the modeling of the 2033 and 1151 cases. 

B. The surface expansion modeling 

In order to model the regressing walls of the 1151 and 2033 channels a surface morphing procedure which 
accounts for the local wall oxidation should be developed. The corresponding procedure described here rests 
on the assumption that initial surface is smooth (no sharp corners) and is defined as a set of adjacent quads. 
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Furthermore, we do not expect the surface mesh to uniquely define the normal direction with respect to a 
particular surface side. The procedure is capable of creating surface meshes of a sufficiently good quality for 
meshing the computational domains and is as follows. 

We start with a control net based channel surface created in Gridpro. We mesh the computational domain 
where the channel surface is one of the boundaries and run the Gasp 5 N-S solver to obtain a CFD solution. 
A stand alone morphing code reads the initial channel surface mesh data and the CFD solution results then 
computes the normal surface vectors as the cross products of the surface element edges. The normals are 
associated and stored with each of the surface element. If nodes are shared among several surface elements 
(which is usually the case, as shown in Fig. 17), an averaging is done to compute the normals associated 
with the surface nodes. Each of the following cross-products OA x OC, OE x OD and OC x OD defines 
a different normal direction at point O. 

The resulting normal direction OF is the result of averaging. We then provide reference points whose 
coordinates are given with respect to the origin of the computational domain. These reference points indicate 
in an absolute sense the inside and outside of the channel. Based on these reference points, the directions 
of the normal vectors associated with each of the surface nodes are adjusted to be uniformly directed with 
respect to the surface side. The problem, that the provided reference points solve, is illustrated in Fig. 18. 
In order to define a surface mesh element it is not sufficient to specify the nodes of the surface with their 
coordinates. In addition to this, the connections between the nodes must also be specified. For example 
let the nodes of a simple surface which consists of only two quads, as shown in Fig. 18, to be counted and 
have the numbers from 1 to 6 as presented in the figure. Then the two quads presented in the figure can, 
for example, be defined as follows: quadl ([1,4] ; [1,2] ; [2,3] ; [3,4]) quad2 ( [4,3] ; [4,6] ; [3,5] ; [6,5] ) . The pairs of 
numbers in the square brackets represent the edges of the quads. Let us now define a universal rule for the 
definition of the normal direction associated with a quad as the cross product of vectors defined by its two 
first edges. This would imply that the normal of the first quad is defined as follows [l,4]x[l,2] and for the 
second quad it will be [4,3]x[4,6]. Note that we assume the edges to be vectors with directions from the first 
nodes of the edges to the second ones. The presented cross-products define the normals to the quads as is 
shown in Fig. 18. Both of the normals are defined correctly but they point to opposite directions. This would 
be a problem in our application of morphing of the surfaces, thus, a special care must be taken to define the 
normals properly. The next step is to find a closest volume node for each surface node with the values of 
the computational results to be associated with the surface nodes. The local oxidation rate at each surface 
node is computed then and stored. With the local oxidation rate and the normal direction obtained, each of 
the surface elements is accordingly moved. The movement of the nodes is done in a time loop incremented 
with the time steps that are much smaller than the time between sequential CFD runs. Visual surface mesh 
quality control achieved by taking the morphed mesh back into the Gridpro and checking that the resulting 
surface is smooth enough. 

If there are imperfections due to the collapsing of the surface elements, the procedure is re-started and 
run up to the time where the surface mesh was still smooth. This intermediate surface mesh serves as a 
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reference for an analytical approximation. After several of experimentations, it was found that a parabolic 
approximation to the channel surface is useful. It is defined as follows. The oxidation rates along the edges 
of the channel (the entrance, the exit and the corners which are defined by the intersection of the channel 
surface with the channel symmetry plane) are computed and stored for each of the rows and columns of the 
channel surface mesh. These values serve to define a series of parabolas associated with each of the rows 
and columns of the surface mesh. Figure 19 shows two examples of parabolic curves for a row and a column 
at the middle of the surface mesh. For this specific case, there are therefore 57 times 238 parabolic surface 
curves which are re-computed every time the surface is morphed, or every 1 sec. 

The surface regression is then performed in two steps. We first morph each surface mesh point by moving 
it half of the local computed CFD defined oxidation rate from the parabolic expression for the row that 
it belongs to. We then repeat the procedure again for each surface mesh point but now move it from the 
parabolic value depending on its associated column. In such a way we take into account the rate along all 
of the channel corners when computing the amount of necessary displacement for each of the surface mesh 
points. This gives a more realistic and a more 3-D meshing friendly approximation of the surface than a 
surface advancement based on the actual local CFD oxidation rates. However this is not the final mesh 
which can be directly used in the 3-D volume meshing, since it is not integrated with the rest of the surfaces 
which define the boundaries of the specimen. The resulting mesh is taken back into Gridpro and is used as 
a control net to produce the spline surface to be integrated with the outer surface of the specimen. Finally 
the spline based surface mesh is used in Gridpro to account for the channel surface mesh change, and a new 
3-D volume computational mesh is created. 

V. Results and Discussion 

A. Case 1151 with delamination meshing and flow results 

Attention was paid to the proper meshing of the computational domain making use of the Gridpro software. 
The multi zone layout of the domain presented in Fig. 20 allowed for the adjustment of the zone boundaries to 
automatically improve the mesh quality in terms of its orthogonality skewness and aspect ratio. To account 
for the possible influence of delamination, a sub-channel was added to the layout of the computational 
domain as shown in Fig. 21. After a preliminary CFD solution, based on a coarser mesh, the position of 
the shock wave was determined and the mesh was adapted to better capture the shock wave gradients as 
presented in Fig. 22. To better resolve the boundary layer gradients mesh points were clustered at the walls 
as presented in Fig. 23. It can be seen that the mesh cells are orthogonal with small aspect ratio (excluding 
the cells in the boundary layer) and with almost no skewness in 80 to 90 % of the cells. 

To resolve mesh singularities at the corners of the specimen and at the entrance of the sub-channel, 
artificial double-sided surfaces were introduced into the computational model at the mesh generation stage. 
After each of the sequential computations, the channel walls and adjacent boundaries (including artificial 
walls) were morphed and the new computational domain was re-meshed in a semi-automated manner. Mesh 
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quality was checked for skewness with tools available in Gridpro for each of the solved cases (see Fig. 24). 
The mesh gradients are extremely large in the cases with a sub-channel since the outer flow, channel flow, 
and the sub-channel flow were all solved in a single run starting from the pre-shock NATA solution. However, 
the mesh quality was never allowed to worsen below the strict limits imposed by the checking procedures. 
After morphing the surface, the re-meshing of the computational domain was relatively straightforward with 
some limited amount of manual work. 

B. Influence of delamination 

In order to understand the behavior of the flow parameters influencing the wall regression rate, temperature 
and atomic oxygen mass fraction are presented in Figs. 25 and 26. Additionally CO mass fraction and velocity 
streamlines are presented in Figs. 27 and 28. Figure 29 presents a comparison between the oxidation rates 
at the channel wall for 0 and 50 seconds. Finally Fig. 30 presents a comparison between the computational 
results for calculations with and without the sub-channel with the experimental data. Figures 25 - 28 show 
that for the expanded sub-channel the flow in the sub-channel entrance area stagnates with recirculation 
areas observed. This results in the atomic oxygen available in the sub-channel area being consumed by the 
chemically active wall. CO, the product of the oxidation, remains in the sub-channel entrance area shielding 
a fresh intake of the atomic oxygen from the channel flow. This in turn results in a lower oxidation rate in 
the sub-channel area compared with that along the main channel wall as can be seen in Fig. 29. Figure 29 
demonstrates that this occurs because the 50 sec surface regression profile is actually reversed compared 
to the initial (0 sec.) result. The phenomena is probably due to two different factors contributing to the 
oxidation rate: pressure build up at the sub-channel corners and the stagnation of the flow in the sub-channel 
area. These factors influence the oxidation rate magnitude in different directions. A more accurate DSMC 
or BGK solution to the sub-channel inlet problem may provide more insight into the phenomena. 

The entrance and exit hole diameters growth, presented in Fig 30, show the influence of the delamination 
(the presence of the sub-channel) on the main channel regression. As can be seen in the figure, the presence of 
the sub-channel increases the rate of the exit diameter growth compared with the results of the computation 
without the sub-channel. This influence brings the computational results closer to the experimental data. 
At the same time, the growth of the entrance diameter of the channel remains the same (does not depend 
on the presence of the sub-channel). 

C. Influence of the flow chemistry upstream of the channel 

Prior to presenting the influence of the flow chemistry upstream of the channel on the wall regression rates 
some explanation must be made about the accuracy of the attempted modeling. Figures 31 and 32 present 
the temperature profiles in the shock wave and along the stagnation line respectively. The locations in 
the flow where the temperature profiles presented in Figs. 31 and 32 are taken from are shown in Figs. 33 
and 34 respectively. As can be seen in the figures the temperature levels at around 9,000 K, a temperature 
sufficiently high that chemical non-equilibrium must be considered in the gas dynamics. 
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Our consideration of the flow chemistry is based on a finite rate model at a single temperature which is 
not very accurate at these temperatures. However we will show that the atomic oxygen generation in the 
shock wave can influence the wall regression rate to some extent. For this reason a comparison of the wall 
regression rates between the N-S solutions with and without chemistry in the flow are presented below. 

Comparisons between the temperatures and mass fractions of atomic oxygen (the two parameters which 
are included into the computation of the wall regression rate) are presented in Figs. 35 and 36 as top 
and bottom figures for cases with and without flow chemistry. Additional comparison of pressures and 
temperatures are presented in Figs. 37 - 38. The set of chemical reactions used in this research is given in 
Table 2 with reaction coefficients being dependent on the temperature as defined in the GASP 5 computational 
tool. The results are provided for an initial channel configuration geometry. It can be seen in the figures 
that the general flow parameters (namely pressure, temperature and x component velocity) have reasonable 
spatial distributions. In particular, flow pressure decreases while the flow accelerates in the channel and 
the temperature of the flow drops because the thermal energy is transformed into the kinetic energy of the 
accelerating flow. It should be noted that some difference exists between the cases with and without flow 
chemistry. The temperature after the shock is smaller in the case with chemical reactions due to dissociation 
reactions which consume thermal energy. The main point with respect to this work is that amount of atomic 
oxygen in the flow is slightly higher in the case with chemical reactions due to the additional production 
of atomic oxygen by dissociation of the oxygen molecules in the shock. Examination of Fig. 26 shows that 
the atomic oxygen mass fraction does not change inside the channel, since there are essentially no chemical 
reactions there, as was discussed earlier in the paper. 

The above phenomena is reflected in Figs. 39 and 40. The somewhat increased amount of the atomic 
oxygen increases the wall oxidation rate which is reflected in the amount of carbon monoxide production as 
presented in Figs. 39 and 40 for the entire computational domain and in the channel area, respectively. A 
comparison between the computational results with and without chemistry and the experiment is presented 
in Fig. 41. For reference, the results for the flow without the sub-channel are also presented. As can be 
seen in the figure, the production of the atomic oxygen in the flow influences the wall regression rate at the 
entrance of the channel, where the flow is exposed to the shock generated atomic oxygen. Unlike the earlier 
result (Fig. 30) which showed that the presence of the sub-channel primarily affects the outlet diameter, here 
the additional atomic oxygen strongly affects the channel entrance region oxidation rate. 

D. Case 2033 meshing and flow results 

Figure. 42 provides a schematics of the three-dimensional wedge case geometry and the layout of the compu- 
tational domain associated with it. Presented in the figure is the 320 second case, with the channel already 
substationally expanded to demonstrate all of the features of the computational domain. The initial channel 
dimensions were however 0.990 x 38.1 mm (0.0387 x 1.5 inches) with a depth of 6 mm. The multi-scale 
external-internal flow for this setup was solved in a single run. The computational domain was structurally 
meshed using the block-structured techniques available in Gridpro with different mesh densities in different 
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blocks in order to accommodate the flow phenomena which ranged from the shock wave structure in front of 
the modeled portion of the Space Shuttle wing leading edge to the boundary layer inside the micro-channel 
crack. The computational domain consisted of the larger outer sub-domain and a tiny well meshed micro- 
channel domain. Figure. 42 also presents the mesh block structure which was optimized to produce a high 
quality 3D mesh covering the entire computational domain. The inlet surface presented in the figure was 
set up with the boundary conditions corresponding to the free stream NATA pre-shock solution presented in 
table 3. The no-slip thermal boundary conditions were set at the surface of the specimen, which is marked 
as wedge in the figure, while the temperature jump velocity slip boundary conditions which accounted for 
the oxidation at the wall were used at the channel walls. First order extrapolation boundary conditions 
were set at the two exits of the computational domain, and two symmetry planes were used to limit the 
computational domain in Z direction. The flow direction is also presented in the figure as a vector in the XY 
plane, making 45 degrees angle with the X axis. Figure 43 presents the computational mesh. As can be see 
in the figure, care was taken to resolve the flow features especially those inside the channel and in vicinity of 
its entrance. Figure 44 further illustrates the features of the computational domain by presenting the Mach 
number contours obtained from a CFD solution with use of the presented computational mesh. 

To better illustrate the channel portion of the computational domain, a mesh block which corresponds 
to the channel and its entrance area is presented in Fig. 45. As can be seen in the figure, several mesh blocks 
were used to mesh the channel area, which required, in part, the boundary layer features resolution, since the 
flow properties at the walls of the channel are essential for the correct prediction of the channel expansion. 
Figure 46 presents the atomic oxygen density contours plotted in three dimensions using the excerpt of the 
computational mesh, presented in Fig. 45. 

To illustrate the 3D shock wave capturing attempted in this research, a portion of the computational 
domain with the exposed interior of the mesh is presented in Fig. 47. This is a 60 second case mesh which 
was used in one of the solutions from the series of solutions covering the expansion of the channel. As can 
be seen in the figure, the mesh resolution is significantly higher in the shock wave area which allows for a 
better shock wave structure modeling. The quality of the solution generated based on this mesh is presented 
in Fig. 48 where the temperature contours are presented in 3D space, showing, in particular, the shock 
wave capturing. The series of N-S calculations was performed based on the computational domains similar 
to the one presented in Fig. 42. The free stream uniform boundary conditions obtained from the NATA 
solution (see Table 3) were applied along the entrance of the computational domain as presented in Fig 42. 
The computational domain accommodates bow and oblique shock waves in front of the specimen. The inlet 
surface is aligned with the shock wave structure to reduce the size of the computational domain. 

A series of iterations at 0, 60, 120, 220, 320 and 420 seconds were performed with the channel shape 
changing according to the oxidation process at the wall. 

Figures 49- 51 show the major flow features in the computational domain, namely Mach number, pressure 
and temperature contours. As can be seen in the figures the overall flow features are well resolved and show 
a reasonable pattern in all portions of the computational domain. In particular, the shock wave structure is 
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resolved and the flow inside the channel is accurately modeled. Figures 52 and 53 show the atomic oxygen 
and CO mass fraction distribution at the channel walls and down stream (over the specimen) for the 60 sec 
case. As can be seen in the figures, a small portion of the product of reaction at the walls of the channel 
escapes the channel through the subsonic boundary layer and moves over the specimen. A corresponding 
decrease in the atomic oxygen mass fraction can also be observed in Fig. 52. Such a mechanism can actually 
reduce the specimen surface destruction in the downstream area since the product of reaction, the CO 
generated inside the channel, when escaping the channel, may shield the the specimen surfaces from the 
chemically active atomic oxygen. 

Figure 54 presents the shape of the channel computed at different times. During the first two computa- 
tional steps, 0 and 60 seconds, the surface morphing procedure was used in full as described in Sec. B. After 
the channel was sufficiently opened with its walls exposed to the coming flow, more reliable values of the local 
oxidation rates become possible to compute. Beginning with 120 case, the parabolic surface approximation 
as defined in section B was omitted and the local surface displacements were computed based on the local 
oxidation rates obtained from the CFD solutions. The rest of the procedure described in section B remained 
the same. 

Following observations can be made from the results of modeling. The initial channel expansion is less 
intensive and more even in all of the directions. However, when the channel opens, which allows more atomic 
oxygen to enter the channel interior, the process accelerates, with the predominant expansion areas being 
formed in the downstream portion of the channel. 

Figure 55 shows a comparison between the channel shapes obtained computationally and experimentally 
in the arcjet facility. The experimentally obtained surface was roughly reconstructed based on the post test 
measurements of the channel. As can be seen in the figure the channel sizes are in a good agreement at 
the entrance of the channel with the experimentally obtained channel being slightly larger. The size of the 
experimental channel at the exit is approximately 20 percent smaller than that of the computational one. 
The overall shapes of the channels are remarkably close. 

VI. Conclusions 

In this work, an approach was developed to model the continuum flow through the carbon crack channel 
of the RCC test samples. The gas dynamics was modeled using the N-S technique with chemically active 
boundary conditions. Atomic oxygen impacting a bare carbon wall is assumed to chemically react with the 
wall to produce CO gas (with a probability of 0.9) that enters the flowfield. It was found that carbon surface 
degradation caused by atomic oxygen accounts for large mass loss and shape similar to the arcjet samples 
for considered cases. 
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Table 1. Free stream flow parameters from NATA arcjet solution and initial channel sizes for case 1151. 


Channel diameter (m) 

0.00325 

Channel length (mm) 

0.006 

Pressure (Pa) 

70.5 

Mach 

5.83 

Temperature (K) 

900 

Chemical species mass fractions 

N 2 = 0.298 


0 2 = 0.2971e-6 


NO = 0.3763e-6 


N = 0.4666 


0 = 0.2345 

Channel wall BC 

Diffuse, oxidation reaction 

Solver 

N-S 
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Table 2. Chemistry. 


Reaction 


Forward rates * 



Equilibrium rates ** 





c 

V 

- e/kT 

c 

V 

- e/kT 

N 2 + N 2 - 

-> N + N + N 2 

4.7 x 10 14 

-0.5 

113000 

18000 

0 

113000 

N 2 + O 2 - 

-4 N + N + 0 2 

1.9 x 10 14 

-0.5 

113000 

18000 

0 

113000 

N 2 + NO - 

-► N + N + NO 

1.9 x 10 14 

-0.5 

113000 

18000 

0 

113000 

N 2 + N^ 

N + N + N 

4.085 x 10 19 

-1.5 

113000 

18000 

0 

113000 

n 2 + o^ 

N + N + O 

1.9 x 10 14 

-0.5 

113000 

18000 

0 

113000 

02 + N 2 - 

■+ 0 + 0 + n 2 

7.2 x 10 15 

-1 

59500 

1,200,000 

-0.5 

59500 

02 + 02 - 

-> 0 + 0 + o 2 

3.24 x 10 15 

-1 

59500 

1,200,000 

-0.5 

59500 

0 2 + NO 

-> 0 + 0 + NO 

3.6 x 10 15 

-1 

59500 

1,200,000 

-0.5 

59500 

o 2 + n^ 

0 + 0 + N 

3.6 x 10 15 

-1 

59500 

1,200,000 

-0.5 

59500 

O 2 + 0 — ■» 

0 + 0 + 0 

9.0 x 10 16 

-1 

59500 

1,200,000 

-0.5 

59500 

NO + N 2 - 

-> N + 0 + N 2 

3.9 x 10 17 

-1.5 

75500 

4000 

0 

75500 

NO + 0 2 

-> N + O + O 2 

3.9 x 10 17 

-1.5 

75500 

4000 

0 

75500 

NO + NO 

^ N + 0 + NO 

7.8 x 10 17 

-1.5 

75500 

4000 

0 

75500 

NO + N- 

> N+O+N 

7.8 x 10 17 

-1.5 

75500 

4000 

0 

75500 

NO + 0 - 

tN+O+O 

7.8 x 10 17 

-1.5 

75500 

4000 

0 

75500 

0 + N 2 -> 

N + NO 

7.0 x 10 4 0 

0 

38000 

4.5 

0 

37500 

NO + 0 - 

■>n + o 2 

3,200,000 

1 

19700 

0.00333 

0.5 

16000 

* Presented 

are the coefficients in 

the formula: Kf 

= CT^eT 

■e/kT 





where: k is Boltzmann constant, T a is the rate controlling temperature defined as follows: 

T a = T“T^T e 1 -“- /3 where: T„ - is vibrational temperature and T e is electron temperature. In the solved cases: a = 1, (3 = 0 
** Presented are the coefficients in the formula: K c = CT 71 e~ e ^ kT 
e/kT is activation energy 
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Table 3. Flow parameters for the case 2033 (obtained from NATA arcjet solution). 


n 2 

0.177 

0 2 

0.033 

N 

0.551 

0 

0.239 

NO 

2.426e-6 

Temperature (K) 

929 

Pressure (Pa) 

151.07 


** Free stream conditions taken from a separate 
arcjet CFD calculation upstream of the arcjet shock. 



Figure 1. Case 1151. The experimental facility. 


NATA 

Analysis 



Arcjet Nozzle 

nozzle size = 10 inches 
Z distance = 8 inches 
mass flow = 0.38 Ibm/sec 
current =1480 amps 
bulk enthalpy = 6529 Btu/lbm 


CFD 
Analysis 
Domain shock 



Figure 2. Case 2033. The schematics of the experimental facility. 
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Figure 3. Initial and post-test specimen configurations. 



Figure 4. The simplest 1151 case topology. 
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Figure 5. The simplest 1151 specimen and domain boundaries. 



Figure 6. 


The simplest 1151 mesh block structure. Lines indicate the morphed topology edges. 
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Figure 7. The simplest 1151 case mesh. Note the red color lines are the topology edges seen in Fig 6. 



Figure 8. The simplest 1151 mesh in the channel Note the red color lines are the topology edges seen in Fig 6 
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Figure 9. The simplest 1151 mesh in the corner at the entrance of the channel. The problems of the mesh 
can be clearly observed. The colored lines except for the grin lines representing the mesh itself represent the 
boundaries of the mesh blocks. 



Figure 10. 1151 improved topology. This is still a 2D topology with all of the new topology nodes lying in the 
same plane with the original topology nodes. 


19 of 41 





Figure 11. 1151 improved topology in the channel region. 



Figure 12. 1151 domain representation with artificial surfaces. 
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Figure 13. 1151 improved block mesh structure. The colored lines represent the boundaries of the mesh 

blocks. The red lines highlight the boundaries of one of the blocks. 



Figure 14. 1151 improved mesh. The entire domain. The colored lines, except for the green lines that represent 
the mesh itself, represent the boundaries of the mesh blocks. 
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Figure 15. 1151 improved mesh. The mesh in the channel 



Figure 16. 1151 improved mesh. The mesh at the corner of the channel entrance. Most of the mesh problems 
are resolved. 
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B C 



Figure 17. The surface normal definition. Point O represents a surface node shared by four surface quadrilat- 
erals, ABCO, CGDO, ODJE and OEHA. Note that point F is not in the plane of any of quads. 




0 



Figure 18. The normal directions with respect to the surface side (with indicated in red) are undefined. 
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Figure 19. The analytical approximation of the channel surface 



-0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08 


X [M] 

Figure 20. Multi zone layout of the computational domain. The presented case corresponds to 50 sec. 
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Figure 21. The domain layout in vicinity of the sub channel. 
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Figure 22. Mesh adaptation in the shock wave area allowed for a better gradient resolution. 


0.005 - 


0.004 



0.002 0.003 0.004 

X [M] 


Figure 23. The sub-channel entrance meshing. Mesh clustering close to the wall is presented along with the 
overall quality of the mesh in vicinity of the sub-channel entrance. The presented case corresponds to 50 sec. 
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Figure 24. Mesh goodness testing information produced by Gridpro. 
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Figure 25. Temperature contours in the channel. [K] 
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Figure 26. Oxygen mass fraction in the channel. 


CO mass fraction 
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Figure 27. CO mass fraction contours in the channel. 
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Figure 28. 


Velocity streamlines. Colors present velocity magnitudes in the channel, [m/s] 



Figure 29. Oxidation rates at 0 and 50 seconds, [m/s] 
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Time, sec 

Figure B-9 Test Data of An Impacted RCC Hole Growth 

Figure 30. Wall regression. Entrance and exit hole diameters. N-S results with and without sub-channel, 
[m/s] 



Figure 31. Temperature in the shock. [K] 
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Figure 32. Temperature along the stagnation line. [K] 
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Figure 33. Location of the shock sample line 
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Figure 34. Location of the stagnation sample line 
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Figure 35. Temperature flow fields with and without flow chemistry [K]. 
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Finite Rate, Velocity Slip at the channel wall, T=Tw 

c+o=co 



Figure 36. Oxygen mass fraction flow fields with and without flow chemistry. 


Finite Rate. Velocity Slip at the channel wall, T=Tw 

c+o=co 



Figure 37. Pressure flow fields with and without flow chemistry [Pa]. 
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Finite Rate, Velocity Slip at the channel wall, T=Tw 

c+o=co 



Figure 38. X component velocity fields with and without flow chemistry [m/s]. 
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c+o=co 



Figure 39. CO mass fraction flow fields with and without flow chemistry. 
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Finite Rate, Velocity Slip at the channel wall, T=Tw 

c+o=co 
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Figure 40. CO mass fraction flow fields in the channel with and without flow chemistry. 



Time, sec 

Figure B-9 Test Data of An Impacted RCC Hole Growth 



Figure 41. Wall regression. Entrance and exit hole diameters. N-S with and without chemistry in addition to 
with and without sub-channel, [m/s] 
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Y 



Figure 42. 320 sec case. Topology of the computational domain. 



Figure 43. 320 sec case. Computational mesh. 
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Figure 44. 320 sec case Mach number contours in the 3D computational domain. 



Figure 45. 320 sec case. An excerpt of the computational mesh in vicinity of the channel area. 
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Density O 



Figure 46. 320 sec case. Atomic oxygen number density in vicinity of the channel. 



Figure 47. 60 sec case. Computational mesh. Illustration of the shock wave capturing. 
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Figure 48. 60 sec case. Temperature contours in the 3D computational domain. 3D shock wave capturing. 



Figure 49. 60 sec case. Mach number contours in the symmetry plane. 
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X[m] 


Figure 50. 60 sec case. Pressure contours in the symmetry plane. 



X[m] 


Figure 51. 60 sec case. Temperature contours in the symmetry plane. 



Figure 52. 60 sec case. O Mass fraction over the specimen. 
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Figure 53. 60 sec case. CO Mass fraction over the specimen. 



Figure 54. Channel shape at different times. 



Figure 55. 


Comparison of the computationally obtained and the the post-test shapes of the channel. 
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